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We study atomistically the fracture of single crystal silicon at atomically sharp notches with 
opening angles of degrees (a crack), 70.53 degrees, 90 degrees and 125.3 degrees. Such notches 
occur in silicon that has been formed by etching into microelectromechanical structures and tend 
to be the initiation sites for failure by fracture of these structures. Analogous to the stress intensity 
factor of traditional linear elastic fracture mechanics which characterizes the stress state in the 
limiting case of a crack, there exists a similar parameter K for the case of the notch. In the 
case of silicon, a brittle material, this characterization appears to be particularly valid. We use 
three interatomic potentials: a modified Stillinger- Weber potential, the Environment-Dependent 
Interatomic Potential (EDIP), and the modified embedded atom method (MEAM). Of these, MEAM 
gives critical isT-values closest to experiment. In particular the EDIP potential leads to unphysical 
ductile failure in most geometries. Because the units of K depend on the notch angle, the shape of 
the K versus angle plot depends on the units used. In particular when an atomic length unit is used 
the plot is almost flat, showing — in principle from macroscopic observations alone — the association 
of an atomic length scale to the fracture process. 
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I. INTRODUCTION 



Recently there has been experimentaliiSi^ and 
theoretical^ interest in fracture in sharply notched sin- 
gle crystal silicon samples. Such samples have techno- 
logical importance because silicon is a commonly used 
material in the fabrication of MEMS devices; the etching 
process used tends to create atomically sharp corners due 
to highly anisotropic etching rates£ Failure in such de- 
vices is often a result of fracture which initiated at sharp 
corners^ In the case of a notch, there exists a parameter 
K analogous to the stress intensity factor of traditional 
fracture mechanics, which parameterizes the elastic fields 
in the vicinity of the notch. Suwito et alm& have carried 
out a series of experiments which have (i) established the 
validity of the stress intensity factor as a fracture crite- 
rion in notched specimens and (ii) measured the critical 
stress intensities for several notch geometries. On the 
theoretical side Zhang 4 has carried out an analysis which 
models the separation of cleavage planes by a simple co- 
hesive law, and thereby derived a formula for the critical 
stress intensity as a function of notch opening angle. The 
material properties which enter this formula are the elas- 
tic constants and the parameters of the cohesive law, the 
peak stress a and the work of separation Tq. This recent 
activity has prompted us to investigate the phenomenon 
of fracture in notched silicon using atomistic simulations: 
In this paper we present direct measurements of the crit- 
ical stress intensity for different geometries (i.e., notch 
opening angles) and compare them to the experimental 
results of Suwito et al. We apply a load by specifying a 
pure if -field of a given strength (stress intensity factor) 
on the boundary of the system. In doing this we are ef- 
fectively using the result of Suwito et al. that the notch 
stress intensity factor is indeed the quantity which deter- 
mines fracture initiation, so we can ignore higher order 
terms in the local stress field. 





FIG. 1: (a) Notch schematic and notation; (b) silicon crystal 
with a notch; darker layer is fixed boundary atoms. 



A. Elastic fields near a notch 

The essential geometry of a notch is shown in Fig. ^ 
The notch opening angle is denoted 7 and the half-angle 
within the material, which is the polar angle describing 
the top flank, is (thus (3 = n — 7/2). As discussed in 
detail by Suwito et al.f2^ it is fairly straightforward to 
solve the equations of anisotropic linear elasticity for a 
notched specimen. The formalism used is known as the 
Stroh formalism, 6 which is useful for dealing with mate- 
rials with arbitrary anisotropy in arbitrary orientations, 
as long as none of the fields depend on the z coordinate 
(this will be the out-of-plane coordinate; note that this 
does not restrict the deformation itself to be in-plane). 
Here we only consider mode I (symmetric) loading. The 
displacement and stress fields for a notch can be written 
as 



Ul = Kr x 9l {e) 



(1) 



(2) 
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FIG. 2: Normal (a) and shear (b) stresses on radial planes as 
functions of plane angle, for 7 = 70°. 



where A plays a role like an eigenvalue; its value is deter- 
mined by applying the traction-free boundary conditions 
to the notch flanks. There is an infinity of possible val- 
ues for A of which we are interested in those in the range 
< A < 1, which give rise to a singular stress field, 
often known as the if-field, at the notch tip. This is 
entirely analogous to the singular field near a crack tip, 
which is simply the limiting case where 7 goes to zero 
((3 — ► 7r), and A becomes one half. Further details of 
the Stroh formalism, as applied to the notch geometry, 
are given in appendix [5] The complete elastic solution 
involves the whole infinity of values for A, corresponding 
to different multipoles of the elastic field. Negative values 
of A correspond to more singular fields which are asso- 
ciated with properties of the core region stemming from 
the non- linear atomistic nature of this region; they do not 
couple to the far-field loading. A > 1 corresponds to fields 
which are less singular, and do not influence conditions 
near the notch-tip, since the displacements and stress 
vanish there. They are, however, essential to represent 
the full elastic field throughout the body, and ensure that 
boundary loads and displacements (whatever they may 
be) are correctly taken into account. This is the basis for 
asserting that only the K- field is important. This field 
is unique among the multipoles in that it both couples 
to the far-field loading and is singular at the notch tip. 
Thus the stress intensity factor must characterize condi- 
tions at the crack tip, and therefore a critical value, K c , 
is associated with the initiation of fracture. The validity 
of this approach hinges on the validity of linear elasticity 
to well within the region in which the K- field dominates. 

From Eq. we see that the units of K and therefore 
K c are stress /length^ -1 which depends continuously on 
the notch angle 7 through A. Hence the shape of a plot 
of K c against notch angle depends on the units used to 
make the plot. In metric/SI units K c changes by an order 
of magnitude between 70° and 125° whereas if an atomic 
scale unit of length is used the plot is nearly flat (Fig.ll5jl. 
The most interesting feature of this is that it seems to 
provide a direct link from macroscopic measurements to 
a microscopic length scale. From a continuum point of 
view, one incorporates atomistic effects into fracture via 
a cohesive zone, a region ahead of the crack tip where 
material cleaves according to a specified force-separation 
law. One of the parameters of such laws is the length 
scale — the distance two surfaces must separate before the 
attractive force goes to zero — which for a brittle material 



TABLE I: Surface energies for silicon according to mSW, 
EDIP and MEAM potentials. 



potential 


surface 


atomic units 


SI units 


mSW 


111 


0.0906 


1.3593 


mSW 


110 


0.1110 


1.6649 


mSW 


100 


0.1570 


2.3545 


EDIP 


111 


0.06538 


1.0475 


EDIP 


110 


0.08194 


1.3128 


EDIP 


100 


0.1320 


2.1150 


MEAM 


111 


0.07668 


1.2285 


MEAM 


110 


0.09030 


1.4469 


MEAM 


100 


0.08126 


1.3019 



is an atomic length scale. It is this scale that one would 
identify from the plot of K c versus angle. Note that one 
can only identify a scale, and not an actual length param- 
eter, in particular because the different geometries that 
are involved in the plot involve different fracture surfaces, 
with presumably different force-separation parameters. 

Fig. [21 shows the normal and shear stresses on radial 
planes (perpendicular to the plane of the sample) ema- 
nating from the notch tip, for unit K and r (i.e., they 
are derived from the tensor /y appearing in Eq. J2J). 
The figures show the functions for the 7 = 70° case; 
the other geometries have the same qualitative behav- 
ior. Both stresses vanish at the maximum angles, corre- 
sponding to the notch flanks; this is in accordance with 
traction- free boundary conditions. What is most impor- 
tant to note is that the normal stress, which presumably 
is most relevant for cleavage on a radial plane, has its 
maximum at 9 = 0. The shear stress, which is relevant 
for possible slip behavior (dislocations) which could com- 
pete with cleavage as a means of relieving stress, is zero 
at 9 = 0, and has a maximum at intermediate angles. If 
there is an easy crystal slip plane in the vicinity of the 
maximum, slip could conceivably compete with cleavage. 

II. SIMULATION 

A. Geometry 

We simulated a cylindrical piece of silicon with a notch, 
making a 'PacMan' shape as in Fig. ^b), consisting of 
an inner core region and an outer boundary region. By 
focusing on just the initiation of fracture we avoid the 
need for large systems since we are not interested in the 
path the crack takes after the fracture (if we were, we 
would have a problem when the crack reached the edge 
of the core region and hit the boundary which is only 
a few lattice spacings away). We consider three notch 
geometries, which we call the 70° (actually 70.5288°), 
90° and 125° (actually 125.264°) geometries respectively, 
referring to the notch opening angles. The 70° sample 
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(a) K = 0.3064 (b) K = 0.3068 

(a) K = 0.3 (b) K = 0.4 



FIG. 3: mSW-crack. 



has {111} surfaces on the notch flanks and the plane 
of the sample is a {110} surface. The 90° sample has 
{110} surfaces on the flanks and the plane is a {100} 
surface (in this case the crystal axes coincide with the 
coordinate axes). The 125° sample has a {111} on the 
bottom flank and a {100} surface on the top flank, while 
the plane of the sample is a {110} plane. In addition, we 
studied the zero degree notch geometry, corresponding 
to a standard crack. The crack plane is a {111} surface 
and is the xz plane in the simulation, and the direction 
of growth is the (211) direction, which is the x direction 
in the simulation. The radius of the inner, core region 
in almost all the cases presented is 5 lattice spacings or 
about 27 A. The exceptions were the crack geometry for 
the EDIP potential (core radius was 7.5 lattice spacings — 
the larger size makes the ductile behavior of the potential 
more obvious) and the 90° geometry with the MEAM 
potential (core radius was four lattice spacings because 
this potential is computationally more demanding). The 
coordinate system in each case is oriented so that the 
plane of the sample is the xy plane and the notch is 
bisected by the xz plane. 



B. Potentials 

We have used three different silicon potentials. The 
first is a modified form of the Stillinger- Weber™ potential 
(mSW) , in which the coefficient of the three body term 
has been multiplied by a factor of 2. This has been noted 
by Hauch et al£ to make the SW potential brittle; they 
were unable to obtain brittle fracture with the unmodi- 
fied SW potential. However it worsens the likeness to real 
silicon in other respects such as melting point and elas- 
tic constantsiSi2iiS The second potential is a more recent 
silicon potential known as "environment-dependent inter- 
atomic potential" (EDIP) ^lii2i which is similar in form to 
SW but has an environmental dependence that makes it 
a many-body potential. Bernstein and coworkers^^^ 
have used EDIP to simulate fracture in silicon. They re- 
ported a fracture toughness about a factor of four too 




(c) K = 0.56 (d) K = 0.6 




(e) K = 0.66 (f) K = 0.76 



FIG. 4: EDIP-crack. 



large when compared with experiment, and that fracture 
proceeds in a very ductile manner, accompanied by sig- 
nificant plastic deformation and disorder. On the other 
hand, using tight-binding molecular dynamics near the 
crack tip they successfully simulated brittle fracture in 
silicon. In view of the failure of many empirical poten- 
tials to simulate brittle fracture, Perez and Gumbsch 16 
used density functional theory to simulate the fracture 
process, measuring lattice trapping barriers for different 
directions of crack growth on different fracture planes. 
A reason for the failure of empirical potentials that has 
been proposed in Ref.ll5lis that their short-ranged nature 
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(a) K = 0.183 



(b) K = 0.1835 
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(a) K = 0.3828 



(b) K = 0.3832 



FIG. 5: MEAM-crack. 



necessarily requires large stresses to separate bonds. This 
however is not the case in our third potential, which is the 
modified embedded atom method (MEAM) of Baskestf 
This is a many-body potential similar to the embedded 
atom method but with angular terms in the electron den- 
sity; it has been fit to many elements including metals 
and semi-conductors. A significant feature of this poten- 
tial is its use of "three-body screening" in addition to 
the usual pair cut-off distance. This means that atoms 
in the bulk see only their nearest neighbors, while surface 
atoms, on the other hand, can see any atoms above the 
surface (for example on the other side of a crack) within 
the pair cut-off distance. The pair cut-off has been set 
to 6 A to allow the crack surfaces to see each other* & 
even after they have separated. The MEAM potential 
has been used successfully to simulate dynamic fracture 
in silicon^ and we have found it to be the most reli- 
able potential in our studies of notch fracture. In table fl] 
we list the low-index (relaxed, unreconstructed) surface 
energies for the three potentials. 



C. Boundary Conditions 

The boundary conditions are as follows: in the z- 
direction (out of the page) there are periodic boundary 
conditions. The thickness of the sample in this direc- 
tion is always one or two repeat distances of the lattice 
in that direction. For the 70° and 125° geometries the 
repeat distance is ^/2a where a is the cubic lattice con- 
stant; for the 90° geometry it is 2a. In the plane, the 
boundary conditions are that an layer of atoms on the 
outside of the system has the positions given by the ana- 
lytic formula Q for displacements from anisotropic lin- 
ear elasticity, with a specified stress intensity factor K. 
The thickness of the layer is twice the cutoff distance of 
the potential, in order that the core atoms feel properly 
surrounded by material^ We interpret the displacement 
formulas in terms of Eulerian coordinates, using an iter- 
ative procedure to compute the current positions. The 
numbers of core atoms were 890, 894, 1260, and 892, for 
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(c) K = 0.4306 



FIG. 6: mSW-70°. 



the 0° (crack), 70°, 90° and 125° systems respectively 
(except for the EDIP/crack case where the core radius 
was 7.5 lattice constants; there the number of core atoms 
was 2002). The number of boundary atoms depends on 
the potential (through the cutoff distance) ; it is typically 
about 500 atoms. For the most part no special consid- 
eration was given to the lattice origin, which meant that 
by default it coincided with the notch tipi^ In a few 
cases it was necessary to shift the position of the origin 
in order to make sure that the notch flanks were made 
cleanly, in particular so that the {111} flanks in the 70° 
and 125° geometries were complete close packed {111} 
surfaces, rather than having dangling atoms. 



D. Critical stress intensities 

The simulation consists of alternating the following 
two steps: (1) We increment the value of if by a small 
amount, changing the positions of the boundary atoms 
accordingly. (2) We relax the interior atoms as follows. 
First we run about 50 steps of Langevin molecular dy- 
namics with a temperature of 500-600 K; the purpose 
of this is to break any symmetry (the 70° and 90° sam- 
ples are symmetric about the xz plane). It is still a zero 
temperature simulation; these finite temperature steps 
are simply a way to introduce some noise. Next we run 
500 time steps of the dynamical minimization technique 



(a) K = 0.2455 



(b) K = 0.2465 



(a) K = 0.266 



(b) K = 0.268 



FIG. 7: EDIP-70 . 



known as "MDmin" (a Verlet time step is carried out, 
but after each velocity update, atoms whose velocities 
have negative dot-products with their forces have their 
velocities set to zero). Finally 500 time steps of conjugate 
gradients minimization are carried out. We observe that 
the combination of both types of minimization is more 
effective (converges to a zero force state more quickly) 
than either alone. The procedure generally results in the 
atoms having forces of around 10 -5 eV/A. 

The initial value for K could be zero; however it turns 
out to be possible to start from a fairly large value of 
K by applying the analytic displacements to the whole 
system at first. When the critical K value, K c , is not yet 
known the increment size is chosen reasonably large to 
quickly find the K c . When this has been found, the sim- 
ulation is restarted from a value below the critical value 
with smaller increments and a more accurate value for 
K c found. The increment is a measure of the uncertainty 
in K c . 

III. RESULTS 
A. Observed fracture behavior 

We observe brittle cleavage of the simulated crystals at 
definite values of K for all geometries using the mSW and 
MEAM potentials, but only for the 70° geometry when 
using the EDIP potential. Figures 151-1141 show snapshots 
of the simulation process for the different geometries and 
potentials. In most cases two or three snapshots are 
shown: one of the configuration immediately before crack 
initiation, one of the configuration immediately after ini- 
tiation, and possibly one of a "late-stage" configuration, 
to illustrate the fracture plane more vividly; generally 
this was chosen to be the configuration corresponding to 
the highest applied load, which depended on how long 
the simulation was run past the initiation point. For the 
EDIP potential, which gives unphysical ductile behavior 
(except in one case, the 70° geometry), more snapshots 
are shown, in order to illustrate the plastic behavior more 



FIG. 8: MEAM- 70°. 



completely, since a variety of stages is involved. 

The behavior in crack geometries is shown in Figs.|3J- 
[5] The initial applied load must be such that no crack 
healing takes place upon relaxation, so that the loca- 
tion of the crack corresponds to the center of the system 
(in reference to which the boundary displacements are 
calculated). In this case we are not investigating crack 
initiation (since the notch is already a crack)but crack 
growth; the critical K c is defined as that at which the 
crack advances, or when the next bond across the crack 
plane breaks. This is somewhat hard to see in the fig- 
ures; one must count atoms along the crack surface and 
compare from one figure to another to see that growth 
has occurred. 

The mSW and MEAM potentials produce similar, brit- 
tle, fracture behavior. The EDIP potential produces 
quite different behavior; the crack propagates in a ductile 
manner. Frame (a) shows the configuration before any 
plastic deformation has taken place. Frame (b) shows 
what appears to be the nucleation of a dislocation onto 
the {110} slip plane which is at an angle of 54.6° to the 
positive x-axis. By frame (c) the crack tip has blunted 
noticeably, and in frame (d) a growth of the blunted crack 
by about a lattice constant has taken place — we take 
the stress intensity at this stage to be the critical value. 
Frames (e) and (f) show a void nucleating and growing 
behind the crack tip, which would under further loading 
join with the crack — coalescence of voids the essence of 
ductile crack growth. 

In the 70° system (Figs. fracture occurs along a 
{111} plane. There are two choices for this, symmetri- 
cally placed with respect to the xz plane. Here all three 
potentials produced brittle behavior; this was the only 
geometry in which the EDIP potential did so. Possible 
reasons for this exception are discussed in section IIVI 
However, when the origin was not shifted as mentioned 
in section 111 Dl so that the notch flanks had dangling 
atoms, the EDIP-behavior was quite different: the notch 
blunted to a width of several atomic spacings. 

The behavior for 90° models is shown in Figs. I9TI11I 
We observe three different behaviors for three different 



(a) K = 0.3860 



(b) K = 0.3863 
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(b) K = 0.491 
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FIG. 9: mSW-90°. 



potentials — providing a cautionary demonstration of the 
limitations of empirical potentials. The easy cleavage 
planes available here are the {110} planes which are ex- 
tensions of the notch flanks. The mSW model starts 
to cleave along the lower of these (the extension of the 
upper flank) but the crack advances only one atomic 
before cleavage switches to an adjacent parallel plane. 
The net result is a kind of "unzipping" along the hard 
{100} plane. This is presumably because the peak in the 
normal stress across this plane, compared to the normal 
stress at the 45° angle, outweighs the increased cost of 
cleavage (but note that the surface energy ratio 7100/7111 
is in fact lower for ME AM, which cleaves on the {110} 
plane — see table |I] for the energies of the different sur- 
faces according to the different potentials). The EDIP 
potential deforms plastically in this depicted in 

the six frames of Fig. ^] It is harder to identify spe- 
cific processes here than in the 70° case, including where 
crack growth starts, though it seems to have definitely 
started by the frame (c)(K c = 0.6). The MEAM poten- 
tial behaves in the manner most consistent with exper- 
iment, namely cleaving on {110} planes, and switching 
from one to the other — this is illustrated dramatically in 
the third frame of Fig. ^] Experimentally, switching be- 
tween planes, when it happens, occurs over longer length 
scales (25/im for the 70° case^), although the behavior at 
atomic length scales has not been examined. Too much 
should not be read into the switching we observe, be- 
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FIG. 10: EDIP-90 . 



cause once cleavage has occurred over such distances the 
proximity of the boundary probably has a large effect on 
the effective driving force on the crack. 

For the 125° geometry (Figs. IT2Tll4[l . there are again 
two {111} planes to choose from but they are not sym- 
metrically placed. Fracture occurs for the mSW and 
MEAM potentials on the one closest to the xz plane, i.e., 
closest to the plane of maximum normal stress, which is 
the (111) plane. The direction of growth is [211], and 
growth proceeds much more readily than in the other 
notch geometries, presumably because it is almost along 
the maximum stress plane. In the EDIP system, plas- 



(a) K = 0.293 



(b) K = 0.2935 



(a) K = 0.43 



(b) K = 0.44 
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FIG. 11: MEAM-90 
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FIG. 12: mSW-125° 
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FIG. 13: EDIP-125°. 



tic deformation is favored over cleavage. This appears to 
proceed as follows: First, slip occurs on the (111) plane 
in the [211] direction, as a single edge dislocation is nu- 
cleated (frame (a)-frame (b)). Next, slip occurs on the 
other {111} plane, the (111) plane, in the [211] direc- 
tion, with two dislocations being nucleated (frame (b)- 
frame (c)-frame (d)), on adjacent (111) planes. In the 
last two frames a void appears and grows. 
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(a) K = 0.2195 (b) K = 0.22 



FIG. 14: MEAM-125°. 



B. Critical stress intensities 

The values of K c , for the different potentials as well 
as from experiment, are listed in table ITT1 and plotted in 
Fig. The increment size for K is listed as an esti- 
mate of the error in K c . The values for ductile fracture 
from the EDIP potential are marked with an asterisk as a 
reminder that the definition of K c in these cases is prob- 
lematic. The experimental value for the crack geometry 
is from Ref. Notice that the critical stress intensi- 
ties for difference angles are almost the same in atomic 
units, and differ by more than a factor of ten in stan- 
dard units^i To check for finite size effects, we repeated 
the measurement on the 70° geometry, but with larger 
radius of 8 A, using the MEAM potential. In this case 
K c was determined to be 0.262 ± 0.001, or about 1.7% 
lower than the value from the smaller system. This indi- 
cates that finite size effects are small, but not negligible. 
To compensate for them without using larger systems a 
flexible boundary method could be used, involving higher 
order "multipoles" of the elastic field, appropriate for the 
notch (i.e., solutions with A < 0), which could be relaxed. 



IV. DISCUSSION 
A. Critical stress intensities 

Comparisons are easier to make when looking at the 
data plotted using atomic scale units. Then the data for 
the two brittle potentials is a gentle, almost horizontal, 
curve. The experimental data mostly lies between that 
for the MEAM potential and that for the mSW poten- 
tial, but significantly closer to the former. The exception 
is the 90° case where the experimental value jumps to 
higher than the mSW value. Since the curves from the 
two potentials are very similar in shape — the main differ- 
ence seems to be an overall shift or factor — and the jump 
in the experimental value at 90° is a departure from this 
shape, it would not be meaningful to assert that the mSW 



TABLE II: Critical stress intensity values for different ge- 
ometries and potentials, including experimental data from 
Refs.EH 



Potential 


Geom. 


K c 


Error 


Griffith 


K c (SI) 


mSW 





0.3068 


0.00036 


0.19509 


4.9 x 10 5 


mSW 


70 


0.3832 


0.00036 


- 


9.6 x 10 s 


mSW 


90 


0.3863 


0.00035 


- 


1.78 x 10 6 


mSW 


125 


0.3304 


0.00016 


- 


1.07 x 10 7 


EDIP 





0.6* 


0.02 


0.14634 


9.6 x 10 5 


EDIP 


70 


0.2465 


0.001 


- 


6.1 x 10 5 


EDIP 


90 


0.5-0.6* 


0.0005 


- 


2.4-2.8 x 10 6 


EDIP 


125 


0.5-0.6* 


0.001 


- 


1.5-1.8 x 10 7 


MEAM 





0.184 


0.0005 


0.16406 


3 x 10 5 


MEAM 


70 


0.2665 


0.0005 




6.57 x 10 5 


MEAM 


90 


0.2935 


0.0005 




1.42 x 10 6 


MEAM 


125 


0.22 


0.0005 




6.47 x 10 6 


Expt 





0.2060 




0.1776 


3.3 x 10 s 


Expt 


70 


0.31 


10% 




7.6 x 10 s 


Expt 


90 


0.43 


10% 




2.1 x 10 6 


Expt 


125 


0.22 


10% 




6.5 x 10 6 



potential does a better job in predicting K c in the 90° 
case. For the other angles the MEAM values are more 
or less within experimental error of experiment: the er- 
ror (standard deviation across all the tested samples) is 
close to 10% in all cases (the error is not available for the 
crack case) , and the percentage differences of the MEAM 
values with respect to the experimental values are —10%, 
-14%, -32% and -0.5% for the 0°, 70°, 90° and 125° 
geometries respectively. The 0.5% is clearly fortuitous. 
Note that the experimental error bar is not enough to 
account for the anomalously high value for the 90° case; 
there must be some feature of the physics or energetics of 
fracture initiation in this geometry that is missing from 
the others, and missing from the simulation. 

An interesting question is why the EDIP potential be- 
haves unlike the other potentials and experiment except 
at one particular geometry, the 70° one. Possibly there 
is some feature of this geometry that suppresses the nu- 
cleation of dislocations. Dislocation Burgers vectors in 
silicon are always in ^(110) directions, since these are the 
shortest perfect lattice vectors in the diamond latticed 
The periodic boundary conditions constrain possible dis- 
location lines to be out of the plane. Moreover, since we 
are considering only mode I and II loading, we expect 
slip to be within the plane, so we consider only edge dis- 
locations. In the 70° geometry the ^(110) direction that 
is available within the plane is at an angle of 90° to the 
x-axis, while the cleavage plane is at an angle of 35.26°. 
Looking at Fig. [3 we can see that the shear stress and 
normal stresses on these planes respectively are both near 
their maximum values, although the ratio of shear stress 
to normal stress is 0.43 (both plots are normalized with 
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FIG. 15: Computed critical stress intensities for the three 
potentials and experiment. 



respect to K and r.). Without knowing the values of 
stress needed to initiate slip or cleavage on these respec- 
tive planes, this ratio is not enough to explain anything. 
What we can do is compare the same ratio for the other 
geometries and check if it is higher in other cases, thereby 
tending to make slip more favorable than cleavage, for 
a given potential (namely, EDIP). The numbers — angles 
for slip and cleavage and the appropriate stress ratio — are 
shown in table ITTll Unfortunately, there is no conclusive 
trend. The ratio for 90° is indeed higher than that for 
70° but the others are lower, and EDIP notches suffer 
plastic deformation in all of the other cases. 

For the crack cases we can make a comparison of our re- 
sults with the so-called Griffith criterion for crack prop- 
agation. This comes from setting the energy release rate 
equal to twice the surface energy. An expression for the 
mode I energy release rate in terms of the stress inten- 



TABLE III: Angles of slip planes and crack planes and ratio 
of shear to normal stress for different geometries. 



geometry 


slip plane 


crack plane 


shear/normal 


70 


90° 


35.26° 


0.43 


90 


45° 


45° 


0.52 


125 


27.34° 


-7.90° 


0.34 





54.6° 


0° 


0.37 



sity factor is given in Ref. [21| setting it equal to twice 
the surface energy leads to the following expression for 
the critical stress intensity factor 



K 



Griffith 



20 



7r&22lm((/!i + ^2)/(MiMa)) 



(3) 



where fi± and [ii are the roots of a characteristic polyno- 
mial which depends on the elastic constants and 622 is an 
element of the compliance tensor for plane strain. The ra- 
tio -ftT c / if Griffith is associated with lattice trapping, when 
fracture is brittle. This ratio is 1.57 for the mSW poten- 
tial and 1.12 for MEAM. These values are respectively 
somewhat larger and somewhat smaller than the ratio 
1.25 determined by Perez and Gumbsch using total en- 
ergy pseudopotential calculations 16 (our K c corresponds 
to their K+). In the EDIP case, where fracture proceeds 
only accompanied by significant plastic deformation, K c 
is four times the Griffith value. 

In our simulations, for a given potential, only one frac- 
ture behavior is observed, in contrast to what was ob- 
served in the experiments of Suwito et al*2> Specifically, 
in the case of the 70° geometry, they observed three dif- 
ferent "modes" (not to be confused with loading modes), 
including propagation on the (110) plane, yet we have 
observed cleavage only on {111} planes in this geome- 
try. It is possible that finite temperature, and the rela- 
tive heights of different lattice trapping barriers, play an 
important role here. More likely it is related to experi- 
mental microcracks or defects near the crack tip. In any 
case, it would be of great benefit to systematically calcu- 
late the barriers for different processes that can occur at 
a notch (or crack) tip, as a function of applied load. 

A further point to note, and a warning, is this: In com- 
paring simulations involving such very small length scales 
(27 A) to experiment it is appropriate to consider the 
question of whether the experimental notches are indeed 
as sharp as we have made our simulated notches. Suwito 
et al^ could only put an upper limit of 0.8/im on the ra- 
dius of curvature of their notches, although notch radii of 
the order of lOnm have been reported in etched silicon. 27 
The addition of just a few atoms right at the notch tip 
would presumably have a significant effect on the ener- 
getics of cleavage initiation. We have not made any in- 
vestigation of this, and this question should be borne in 
mind given the absence of experimental data character- 
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izing the notch tip at the atomic scale. Nevertheless, the 
success of our simulations provides an important indica- 
tion that these notches are indeed atomistically sharp. 

V. SUMMARY 

We have determined by atomistic simulation the criti- 
cal stress intensities to initiate fracture in notched single 
crystal silicon samples. The samples had angles of 0° 
(a crack), 70.5233°, 90° and 125.264°— chosen so that 
the flanks of the notches were low index crystal planes. 
These geometries correspond to those studied experimen- 
tally in measurements of critical stress intensities for frac- 
ture initiation. Of the three potentials used, modified 
Stillinger- Weber (mSW), environment-dependent inter- 
atomic potential (EDIP) and modified embedded atom 
method (MEAM), MEAM produced the most realistic 
behavior. The mSW potential produced brittle fracture, 
but its resemblance to silicon in other respects is quite 
weak. Except in the case of the 70° notch, the EDIP 
potential gives ductile fracture with a critical stress in- 
tensity factor, which is much higher than that determined 
using the other potentials, or by experiment. 
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and the length unit to be a = 2.0951 A, the authors 
modeled molten silicon^ However other authors22 have 
taken the energy unit to be e = 2.315eV. The differ- 
ence is not really important since we have modified the 
potential itself to make it more brittle so the resem- 
blance to real silicon is reduced noticeably. When ex- 
pressing quantities in terms of eV— Aunits we use the 
second scaling which is more common. The EDIP and 
MEAM potentials have e = leV and a — lAbuilt in 
as their units. Since a ~ Kr x ~ 1 , the units of K are 
[stress] /[length]^ 1 = [energy]/[length] 2+x , so to con- 
vert a value for K in atomic units to SI units, one uses 
the conversion factor e/a 2+x . Table llVl gives the factors 
for the three potentials and the geometries studied in this 
paper. 



TABLE IV: Unit conversion factors for K. 



Potential 


geometry 


A 


factor 


mSW 





0.5 


1602000 


mSW 


70 


0.51954 


2510000 


mSW 


90 


0.54597 


4620000 


mSW 


125 


0.63047 


32320000 


EDIP 





0.5 


1602000 


EDIP 


70 


0.51922 


2490000 


EDIP 


90 


0.54708 


4730000 


EDIP 


125 


0.62844 


30840000 


MEAM 





0.5 


1602000 


MEAM 


70 


0.51875 


2467000 


MEAM 


90 


0.54794 


4832000 


MEAM 


125 


0.62639 


29420000 



APPENDIX A: UNITS AND CONVERSIONS 

Three different sets of units are used in this paper. To 
each atomic potential (Stillinger- Weber, EDIP, MEAM) 
is associated a set of atomic units (EDIP and MEAM 
use the same units); also we often wish to use SI units to 
compare to experiment. In the context of this paper there 
is the further subtlety that the units of the chief quantity 
under consideration, namely the stress intensity factor K, 
are not simple powers of base units but involve a non- 
trivial exponent A which is a function of geometry and 
potential. In fact the SI units units for K are Pam 1_A 
which for brevity we simply refer to as 'standard units' 
in the paper. 

The units for an atomic potential are determined by 
specifying the unit of energy and that of length (for dy- 
namics the unit of time is determined from these and 
the particle mass). The SW potential as originally writ- 
ten down did not have units built into it. By taking 
the energy unit to be e = 2.1672eV = 3.4723 x 10~ 19 J 



APPENDIX B: STROH FORMALISM FOR 
NOTCHES 

Here we summarize the application of the Stroh for- 
malism to th e notch problem. More details are available 
in Refs. 13*3. 23. We can write the solution for the dis- 
placement field u and the stress function <fr as 

6 

u = ^a Q / Q (z Q ) (Bl) 

a=l 



6 

<j>=^b a f a {z a ) (B2) 

a=l 

The independent variable here is the complex variable 
z a = x\ + p a X2- The stress function <fi determines the 
stresses through an — —4>i,2 and Oii = The p a ,a. a 
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and b Q come from solving the following eigenvalue prob- 
lem: 



satisfied on the bottom flank 9 = —(3. Applying the 
condition on the top flank leads to a matrix equation 



(Q+p(R + R T )+]9 2 T)a = 



(B3) 



where 





C n 


Ci6 


C15 


Q = 




C66 


C56 




Cl5 


C56 


C55 




Cm 




C46 


T = 




C22 


C24 




C46 


C24 


C44 



R 



Ci6 C12 C14 
C66 C26 C*46 
C56 C*25 C45 



(B4) 



The above is general within the context of two- 
dimensional anisotropic elasticity. To specify the notch 
problem we choose a form of the arbitrary function / 
to which we can apply the boundary conditions of the 
problem — that notch flanks are traction-free. The fol- 
lowing choice does the trick: 



b Q q 



(B5) 

where £(0) = cos(9)+p a sin(0) and q is to be determined. 
The traction with respect to a radial plane at angle 9 is 
given by 



t = £ 



b Q b T Q q = -(, 

r 



(B6) 



With the above form the traction condition is already 



K(A)q(A) = 



(B7) 



The appropriate value of A is determined by setting the 
determinant of the matrix equal to zero and solving the 
resulting equation numerically. In the range < A < 1, 
two values can be found, corresponding to modes I and 
II, A 7 and A 77 . For a given A, the vector q is determined 
up to a normalization which is related to how one defines 
the stress intensity factor K. Thus we obtain expressions 
for the displacements which are used in the simulation to 
place the boundary atoms. In the crack case, because X 1 
and A 77 are degenerate at the value 1/2, the definition of 
modes I and II is a little subtle. The {111} plane is not 
a plane of symmetry of the cube, and thus one cannot 
expect to separate modes by their symmetry properties 
as in for example, the isotropic case; following Ref. 0, 
mode I is defined as that for which 1712(0 = 0) = and 
mode II that for which (722(6 = 0) = 0. 

For the purpose of the simulations described in this 
paper, we calculated the Stroh parameters as follows. 
For each potential, the elastic constants were determined 
by standard methods (straining the supercell, relaxing, 
measuring the relaxed energy per unit undeformed vol- 
ume and fitting to a parabola). This gives Cn, C12, and 
C44, which are the three independent constants for a cu- 
bic crystal. In the formulas for the displacements and 
stresses given above, the coordinate system is aligned 
with the notch (in that the negative x-axis bisects the 
notch itself) and not with the crystal axes. So we must 
transform the elastic constants accordingly. Once we 
have the transformed constants we can construct the 
Stroh matrices Q, R and T, and compute the Stroh 
eigenvalues and eigenvectors as above. 
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When applying singular elastic deformations, a check en- 
sures that an atom sitting at the location of the singularity 
is simply not displaced. 

Note that the exact conversion factor depends on the eigen- 
value A which depends on the potential (see appendix |A"1. 
but for a given angle the dependence on potential is quite 
small. 

This fact is mentioned without any citation in Ref. 0. 



